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Abstract 

This paper introduces an extension of the time-splitting sine-spectral (TSSP) 
method for solving damped focusing nonlinear Schrodinger equations (NLS). 
The method is explicit, unconditionally stable and time transversal invariant. 
Moreover, it preserves the exact decay rate for the normalization of the wave 
function if linear damping terms are added to the NLS. Extensive numerical 
tests are presented for cubic focusing nonlinear Schrodinger equations in 2d 
with a linear, cubic or a quintic damping term. Our numerical results show 
that quintic or cubic damping always arrests blowup, while linear damping 
can arrest blowup only when the damping parameter 5 is larger than a thresh- 
old value (5th- We note that our method can also be applied to solve the 3d 
Gross-Pitaevskii equation with a quintic damping term to model the dynamics 
of a collapsing and exploding Bose-Einstein condensate (BEC). 
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1 Introduction 



Since the first experimental realization of Bose-Einstein condensation (BEC) in di- 
lute weakly interacting gases the nonlinear Schrodinger equation (NLS) has been 
used extensively to describe the single particle properties of BECs. The results ob- 
tained by solving the NLS showed excellent agreement with most of the experiments 
(for a review see [4, 12, 11]). In fact, up to now there have been very few exper- 
iments in ultracold dilute bosonic gases which could not be described properly by 
using theoretical methods based on the NLS [22, 25]. 

Recent experiments by Donley et al. [13] provide new experimental results for 
checking the validity of describing a BEC by using the NLS in the case of attractive 
interactions (focusing nonlinearity) in 3d. Since the particle density might become 
very large in the case of attractive interactions inelastic collisions become important 
and cannot be neglected. These inelastic collisions are assumed to be accounted for 
by adding damping terms to the NLS. Two particle inelastic processes are taken 
into account by a cubic damping term while three particle inelastic collisions are 
described by a quintic damping term. Collisions with the background gas and feeding 
of the condensate can be studied by adding linear damping terms. One of the 
major theoretical challenges in comparing results obtained in the experiment with 
theoretical results is to find reliable methods for solving the NLS with a focusing 
nonlinearity and damping terms in the parameter regime where the experiments are 
performed. 

The aim of this paper is to extend the time-splitting sine-spectral method (TSSP) 
for solving the focusing NLS with additional damping terms and to present extensive 
numerical tests. The comparison of our numerical results with the experimental 
results obtained for a collapsing BEC [13] will be presented elsewhere [9]. 

We consider the NLS [8, 38] 

iil>t = -- Aip + V(x) tp - f3\ip\ 2a ip, t>0, x G M. d , (1.1) 
^(x, t = 0) = ^o(x), x G R d , (1.2) 

with a > a positive constant, where a — 1 corresponds to a cubic nonlinear- 
ity and cr = 2 corresponds to a quintic nonlinearity, V(x) is a real-valued po- 
tential whose shape is determined by the type of system under investigation, and 
(3 positive/negative corresponds to the focusing/defocusing NLS. In BEC, where 
(1.1) is also known as the Gross-Pitaevskii equation (GPE) [35], ip is the macro- 
scopic wave function of the condensate, t is time, x is the spatial coordinate and 
V(x) is a trapping potential which usually is harmonic and can thus be written as 

V(x) = \ {^i\x\ + V 7^1) with 71, • • • , 7d > 0. Two important invariants of (1.1) 

are the normalization of the wave function 

N(t) = [ |V(x, t)| 2 dx, t>0 (1.3) 

JR d 
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and the energy 



E(t) 



\ | W(x, t) | 2 + V(x) |^(x, t)| 2 ^— |^(x, t) | 2ct + 2 

Z (J + 1 



dx, t > 0. 



(1.4) 

From the theory for the local existence of the solution of (1.1), it is well known 
that if t)\\m is bounded, the solution exists for all t [38]. As a result, when the 
NLS is defocusing ((3 < 0), conservation of energy implies that J Rd | V^(x, t)\ 2 dx is 
bounded and the solution exists globally. On the other hand, if the NLS is focusing 
(f3 > 0) at critical (ad = 2) or supercritical (ad > 2) dimensions and for an initial 
energy E(0) < 0, the solutions of (1.1) can self-focus and become singular in finite 
time, i.e. there exists a time t* < oo such that [38] 

lim |Vt/>|z,2 = oo and lim = oo. 

However, the physical quantities modeled by ip do not become infinite which im- 
plies that the validity of (1.1) breaks down near the singularity. Additional physical 
mechanisms, which were initially small, become important near the singular point 
and prevent the formation of the singularity. In BEC the particle density \ip\ 2 be- 
comes large close to the critical point and inelastic collisions between particles which 
are negligible for small densities become important. Therefore a small damping (ab- 
sorption) term is introduced into the NLS (1.1) which describes inelastic processes. 
We are interested in the cases where these damping mechanisms are important and, 
therefore, restrict ourselves to the case of focusing nonlinearities (5 > 0, where j3 may 
also be time dependent. We consider the following damped nonlinear Schrodinger 
equation: 

iA = ~ A^ + \/(x)^-/5|^| 2 >-^(|V| 2 )^, t>0, x G M. d , (1.5) 
^(x, t = 0) = ^o(x), x G M d , (1.6) 

where g(p) > for p = \ip\ 2 > is a real- valued monotonically increasing function. 

The general form of (1.5) covers many damped NLS arising in various different 
applications. In BEC, for example, when g(p) = 0, (1.5) reduces to the usual GPE 
(1.1); a linear damping term g(p) = S with 5 > describes inelastic collisions with 
the background gas; cubic damping g(p) = 5\j3p with 5± > corresponds to two- 
body loss [37, 36]; and a quintic damping term of the form g(p) = 5 2 (3 2 p 2 with S 2 > 
adds three-body loss to the GPE (1.1) [37, 36]. It's easy to see that the decay of 
the normalization according to (1.5) due to damping is given by 

N '^ = 7l / hKx,t)| 2 dx = -2/ ^(|^(x,t)| 2 )|^(x,t)| 2 dx<0, t>0. (1.7) 
at JR d JR d 

Particularly, if g(p) = 5 with 8 > 0, the normalization is given by 

N(t)=f \^(x,t)\ 2 dx = e~ 26t N(0) = e - 2St [ |^ (x)| 2 dx, i>0. (1.8) 

JR d JR d 
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There has been a series of recent studies which deals with the analysis and 
numerical solution of the damped NLS. Fibich [16] analyzed the effect of linear 
damping (absorption) on the critical self-focusing NLS, Tsutsumi [39, 40] studied 
the global solutions of the NLS with linear damping, the regularity of attractors 
and approximate inertial manifolds for a weakly damped NLS were given in Goubet 
[21, 20] and by Jolly et al. [26]. For numerically solving the linearly damped NLS 
Peranish [34] proposed a finite difference scheme and this method was revisited 
recently by Ciegis et al. [10] and Zhang et al. [41]. Moebs [32] presented a multilevel 
method for weakly damped NLS and applied it to solve a stochastic weakly damped 
NLS in [31]. Variable mesh difference schemes for the NLS with a linear damping 
term were used by Iyengar et al. [24]. 

Also the TSSP, which we will use in this paper to solve the damped NLS, was 
already successfully used for solving the Schrodinger equation in the semiclassi- 
cal regime and for describing Bose-Einstein condensation using the Gross-Pitaeskii 
equation by Bao et al. [5, 6, 8]. The TSSP is explicit, unconditionally stable and 
time transversal invariant. Moreover, it gives the exact decay rate of the normaliza- 
tion when linear damping is applied to the NLS (i.e. g(p) = 5 with 5 > in (1.5)) 
and yields spectral accuracy for spatial derivatives and second-order accuracy for the 
time derivative. Thus this method is a very good candidate for solving the damped 
NLS, especially in 2d or 3d. We test the novel numerical method extensively in 2d. 

Finally, we want to emphasize that the NLS is also used in nonlinear optics, 
e.g., to describe the propagation of an intense laser beam through a medium with 
a Kerr nonlinearity [18, 38]. In nonlinear optics ip = ^(x, t) describes the electrical 
field amplitude, t is the spatial coordinate in the direction of propagation, x = 
(xi, ■ ■ ■ , Xd) T is the transverse spatial coordinate and V(x) is determined by the 
index of refraction. Nonlinear damping terms of the form g(p) = 5f3 q p q with 5, q > 
correspond to multiphoton absorption processes [16]. 

The paper is organized as follows. In section 2 we present the time-splitting sine- 
spectral approximation for the damped nonlinear Schrodinger equation. In section 3 
numerical tests are presented for the cubic focusing nonlinear Schrodinger equation 
in 2d with a linear, cubic or quintic damping term. In section 4 some conclusions 
are drawn. 

2 Time-splitting sine-spectral method 

In this section we present a time-splitting sine-spectral (TSSP) method for solv- 
ing the problem (1.5), (1.6) with homogeneous periodic boundary conditions. For 
simplicity of notation we shall introduce the method for the case of one spatial di- 
mension (d = 1). Generalizations to d > 1 are straightforward for tensor product 
grids and the results remain valid without modifications. For d — 1, the problem 
becomes 

iipt = ~4> xx + V(x)ij-f3\ij\ 2 ^-igm 2 )ij, a<x<b, t > 0, (2.1) 
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il>(x, t = 0) = ip (x), a<x <b, i>(a, t) = ip(b, t) = 0, t > 0. (2.2) 



2.1 General damping term 

We choose the spatial mesh size h = Ax > with h — (b — a)/M and M an even 
positive integer, the time step is given by k = At > and define grid points and 
time steps by 

Xj := a + j h, t n := n k, j = 0, 1, • • • , M, n — 0, 1, 2, • • • 

Let ip™ be the numerical approximation of ip(xj,t n ) and ip n the solution vector at 
time t — t n — nk with components ■0". 

From time t — t n to time t = t n +i, the damped nonlinear Schrodinger equation 
(2.1) is solved in two steps. One solves 



i i/, t (x, t) = V(x)ij(x, t) - P\ij(x, t) | 2 ^(x, t) - % g(\ij(x, t)\ 2 )ij(x, t), (2.4) 



again for the same time step. Equation (2.3) is discretized in space by the sine- 
spectral method and integrated in time exactly. For t e [t n ,t n+ i], multiplying the 
ODE (2.4) by ip(x,t), the conjugate of ip(x,t), one obtains 



= V(a:)|^(a:^)| 2 -/3|^(a;,t)| 2CT+2 -^(|^(a;,t)| 2 )|^(a;,i)| 2 . (2.5) 



Subtracting the conjugate of Eq. (2.5) from Eq. (2.5) and multiplying by —i one 



-\ij(x,t)\ 2 = fo(x,t)il)(x,t)+ii> t (x,t)ii>(x,t) = -2g{\ij{x,t)\ 2 )\i>{x,t)\ 2 . (2.6) 




XXI 



(2.3) 



for one time step, followed by solving 



obtains 



Let 




f- 1 (f(s) 
0, 



2r) , s > 0, r > 0, 
s = 0, r > 0. 



(2.7) 



Then, if g(s) > for s > 0, we find 



< h(s, t) < s, 



for s > 0, r > 



(2.8) 



and the solution of the ODE (2.6) can be expressed as (with r = t — t n ) 



0<p(t) = p(t n + r):=\^(x,t)\ 2 = h{\iP(x,t n )\ 2 ,t-t n ):=h(p(t n ),T) 
< p(t n ) = |^(a;,t n )| 2 , t n <t< t n+1 . 



(2.9) 
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Combining Eq. (2.9) and Eq. (2.4) we obtain 

i^tOM) = V(x)^(x,t)-p[h(j^(x,t n )\ 2 ,t-t n )Y^(x,t) 

-ig(h(\iP(x,t n )\ 2 ,t-t n ))^(x,t), t n <t<t n+1 . (2.10) 
Integrating (2.10) from t n to t, we find 

tP(x, t) = exp {i [-V(x)(t -t n ) + G (\^(x, t n ) | 2 , t - t n )] - F (\^(x, t n )\ 2 , t - t n )} 
x^(x,t„), t n <t<t n+1 , (2.11) 

where we have defined 

F(s,r) = r 9 (h(s,T)) dr > 0, G(s,r) = T (3 [h(s : r)] a dr. (2.12) 
Jo Jo 

To find the time evolution between t — t n and t = t n+ i, we combine the splitting 
steps via the standard second-order Strang splitting (TSSP) for solving the damped 
nonlinear Schrodinger equation (2.1). In detail, the steps for obtaining from 
are given by 

= exp {-F (\^\ 2 , k/2) + i [-V( Xj )k/2 + G (\^\ 2 , k/2)} } ^, 

M-l 

^** = j2 e -^?/2 ft sin ^ Xj _ a )) 5 j = 1, 2, • • • , M — 1, (2.13) 



= exp {-F k/2) + i [-V( Xj )k/2 + G (|^T, k/2)] } Vf, 

where Ui are the sine-transform coefficients of a complex vector U = (U , U 1: • • ■ , U M ) 
with Uq = Um = which are defined as 



Af-l 



7r/ - 2 

^ = &TT^' ^ = m ^ sin (M^ - a)), ^ — 1) 2, • • • , M — 1, 



where 



$ = il>(x j ,0) = Mxj), J = 0,1,2, 



M. 



(2.14) 



(2.15) 



Note that the only time discretization error of TSSP is the splitting error, which is 
second order in k if the integrals in (2.7) and (2.12) can be evaluated analytically. 



2.2 Most frequently used damping terms 

In this subsection we present explicit formulae for using TSSP when solving the NLS 
with those damping terms most frequently appearing in BEC and nonlinear optics. 

Case I NLS with a linear damping term. We choose g(p) = 5 with 5 > in (1.5). In 
BEC this damping terms describes inelastic collisions of condensate particles with 
the background gas. From (2.7), we find 

f(s) = [ —ds = -Ins and h(s, r) = e - 2Sr s. (2.16) 
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Substituting (2.16) into (2.9) and (2.12), we obtain 

p{t) = e- 2 ^-*") \if>(x, t n )\\ t n <t< t n+1 , 
F(s, r) = 5r, 

G(a>r) = g(l- C — ). 



(2.17) 
(2.18) 

(2.19) 



Substituting (2.18) and (2.19) into (2.13), we get the following second-order time- 
splitting sine-spectral steps for the NLS with a linear damping term 



= exp {-k6/2 + i [-V( Xj )k/2 + (3\^\ 2a (l - e~ 5ak ) /(25a)] } ^, 



M-l 



3 



]T e-^/ 2 ^* sm^xj-a)), j = 1, 2, • • • , M - 1, 



(2.20) 



n+l 



exp 



{-Arf/2 + % [-V( Xj )k/2 + p\iJ**\ 2rT (l - e' 5ak ) /(25a)] } ^ 



Case II NLS with a damping term of the form g(p) = 5[3 q p q , where 5, q > in 
(1.5). For q = 1 (q = 2) we obtain the damping term describing two (three) particle 
inelastic collisions in BEC. From (2.7) we get 



/(S) = /^ 



-ds 



and h(s, r) 



d3"»i- r ~ qSpisi "' (1 + 2q5T(3 q s q ) 1/q ' 

Substituting (2.21) into (2.9) and (2.12), we obtain 

u ix,t n )\ 2 



(2.21) 



P{t) [l + 2q5(3*(t-t n )\i,(x,t n )H lq ' 

F(s,r) = ^-ln(l + 2q5r(3 q s q ), 
2q 

( gi-q 

In (1 + 2q5rf3 q s q ) 



t n <t<t 



n+l, 



G(s,r) = < 



25q 

pi-q s *-q _l + (l + 2q5rp q s q Y q - a 



)/<! 



25(q - a) 



q = a, 



a^q. 



(2.22) 
(2.23) 

(2.24) 



Substituting (2.23) and (2.24) into (2.13), we get the following second-order time- 
splitting sine-spectral method for the NLS 



' exp {i [-V(xj)k/2 + I3 x -i In (l + 6qk/3 c '\^\ 2c ') /(2Sq)] } 



V = 



exp < i 



(I + q5kf3<i\^\ 2 i) 



l/2q 



+ 



2&{q-a) 



-1+ (l + Sqk^l^l 2 ")' 



(1 + qSkpi\^\ 2 i) 



l/2q 



a = q, 



M-l 



*r = E 



ikfj,f /2 



■0;* sin(fj,i(xj - a)), j = 1, 2, • • • , M - 1, 



(2.25) 
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f exp {i [-V( Xj )k/2 + In (l + 5qkp>\il>?\ **) /(2Sg)] } 



n+l 



exp 



{• 



(i + ^fc/39|^;*| 2 9) 



l-?|,/,**|2o— 2 ? 



1/2? 



+ p 2 ^-' CT) - 1 + ( x + *«wri 2 ')" 



(1 + qSk^**^ 



1/2? 



a = q, 



V>**> ° Pl- 



ease III Focusing cubic NLS with a damping term that accounts for two-body and 
three-body loss in a BEC [37], i.e., we choose a — 1, g(p) — 5±(3p + bifi 1 p 2 with 
6i, 5 2 > 0, in (1.5). Using (2.7), we get 



.^L_ + |ln(5 2 /3 + 5 1 /s), s>0, 



(2.26) 



[0, s = 0. 

Substituting (2.7) into (2.12) and changing the variable of integration we obtain 



F(s,r) = J o r g(r 1 (f(s)-2r)) dr** 



=(f(s)-f(h))/2 rKs,r) 1 



£' ~g{h)f{h)dh 



fh{s,r) 



dh 



-\\n{h(s,r)/s) , s > 0, 



s = 0; 



/ 2/i 
where /i(s, r) is the solution of 

/(«) - f(H s , r )) = 2r > for an y r > °> 

with / given in (2.26). Similarly we find 

G(s,r) = J 



dh — I 2Sl s(S 1 +S 2 i3h(s,r))> 

2g(h) \ 0, s = 0. 



(2.27) 



(2.28) 



(2.29) 



Substituting (2.27) and (2.29) into (2.13) we get the following second-order time- 
splitting sine-spectral steps for the NLS with a combination of cubic and quintic 
damping terms 



i^i 



■ exp < I 



V( Xj )k 1 ^h(\^\ 2 ,k/2)(S 1 +S 2 (3\^\ 2 ) 



25i |^| 2 (<5i+«5 2 /3/ l (|V™| 2 ,fc/2)) 



I 0, 

M-l 



= o, 



= J2 e ^' /2 ^i sin (wfe-")), j = l,2,---,M-l, 



■ N /MiV'ri 2 .fc/2) 



exp < % 



V( Xj )k 1 , MIV>**| 2 ,fc/2)(<5i +<5 2 /3|V**| 2 ) 



■ In 



2<5i |^;*| a («i+* a /9h(l^"l a .*/2)) 



(2.30) 



V>" = o. 
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Remark 2.1 As demonstrated in this subsection, the integrals in (2.7) and (2.12) 
can be evaluated analytically for the damping terms which most frequently appear 
in physical applications. If the integrals in (2.7) or (2.12) can not be evaluated 
analytically or the inverse of f in (2.7) can not be expressed explicitly, e.g., if g{p) 
in (1.5) is not a polynomial, one can solve the following ODE numerically by either 
second- or fourth- order Runge-Kutta method 



to get h(s, k/2) for any given s > and set h(s, k/2) = for s = 0. By changing 
the variable of integration in (2.12), see detail in (2.27) and (2.29), the first integral 
in (2.12), i.e. F(s,k/2), can be evaluated exactly (see detail in (2.27)), and the 



numerically by using a numerical quadrature, e.g., the trapezoidal rule or Simpson's 
rule. 

The scheme TSSP is explicit and is unconditionally stable as we will demonstrate 
in the next subsection. Another main advantage of the time-splitting method is its 
time transversal invariance, which also holds for the NLS and the damped NLS 
themselves. If a constant a is added to the potential V, then the discrete wave 
functions ipj' n+1 obtained from TSSP get multiplied by the phase factor e -' ia ( n + 1 ) k ^ 
which leaves the discrete normalization unchanged. This property does not hold for 
finite difference schemes. 

Remark 2.2 For the focusing cubic NLS with a quintic damping term describing 
three-body recombination loss and an additional feeding term for the BEC [27] we 
choose a —\, g{p) = —Si + 52j3 2 p 2 with 61,62 > in (1.5). The idea of construct- 
ing the TSSP is also applicable to this case although we could not prove that it is 
unconditonally stable due to the feeding term. Inserting the above feeding term into 
(2.7), we get 



dhjt) 
dt 



2g(h(t)) h(t), 



0<t< k/2, 



h(0) = s, 




dh, can be evaluated 




(2.31) 



Inserting (2.31) into (2.9), we find 




(2.32) 



and substituting (2.32) into (2.9) and (2.12), we obtain 

mx,t n )\ 2 V6~i 




(2.33) 
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F(s, r) = -Sir + - In [l + S 2 f3 2 s 2 {e^ r - 1)/^ 



1 ps^5~ 2 e 2rh + Js~i + S 2 f3 2 s 2 (e 4 ^ - 1) 
G(s,r) = - — ;===ln 



(2.34) 
(2.35) 



Inserting (2.34) an d (2.35) into (2.13), we get the following second-order time- 
splitting sine-spectral steps for the NLS with a quintic damping term and a feeding 
term 



= fe<5i/2 



cxp 



In 



P\^\ 2 ^e ks i+ y JS l +S 2 p\^\*{e^i-l) 



I 1/4 



[l + £ 2 /? 2 |^| 4 (e 2fe ^ - l)/tfi] - 

M-l 

= E e ~' fcAl ' 72 ^ sin («fe - «)), i = 1, 2, • • • , M - 1, 

z=i 



(2.36) 



e fe5!/2 exp 



n+l 



P\i>" I 2 Vfee" 1 + ^<5l+<5 2 /3 2 |^**| 4 (e 2 ' c ' 5 i -l) 



V(x 3 )fe l i 

' ! 2 h 2V<5l< s 2 m VSl+P\4>"\ 2 Vh 



[l + £ 2 /3 2 |V** | 4 (e 2fe5 i -lj/tfi] 



1/4 



Remark 2.3 T/ie scheme TSSP (2.13) can easily be extended for solving the com- 
plex Ginzburg- Landau equation (CGL) [17, 30] 



iPt = - (1 - i e) - H=V - i (W| 2 - 5i) V>, 



(2.37) 



where e, Si and 5 2 are positive constants. The idea of construncting the TSSP for 
the damped NLS is also applicaple to the CGL provided that we solve 



i tp t = - (1 - i e) Aip, 



(2.38) 



in the first step instead of (2.3). Inserting a — 1, (5 — 1 and g(p) = 5 2 p — Si with 
Si, S 2 > into (1.5) and using (2.7) we get 



m = | 

Inserting (2.39) into (2.7) we find 

h(s,r) 



£ln|<S 2 -<$i/s|, s>0, 
0, s = 0. 



s5i 



S5 2 (l -e- 2 ^) +(J ie -2r«i' 

and substituting (2.40) into (2.9) and (2.12) we obtain 

Si \^{x,t n )\ 2 



tit) 



S 2 \ijj(x,t n )\ 2 (I - e~ 2 ^) + Sie~ 2 ^ 



, t n ^ t ^ ^n+i, 



F(s, r) — —-la 



Si 



2 s<5 2 + {Si - s5 2 ) e- 2 ^ ' 
nl , 1 , 5i - sS 2 + s<5 2 e 2 ^ 

G(s ' r) = 2^ ln si • 



(2.39) 

(2.40) 

(2.41) 
(2.42) 
(2.43) 
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Inserting (2.42) and (2.43) into (2.13), we get the following second-order time- 
splitting sine-spectral steps for the CGL (2.37) 



Si 



exp 



— -In 

2S 2 



Si - S 2 \4>?\ 2 + 6 2 \i%\ 2 e kS * 
5i 



M-l 



= E e- (£+<)fe ^ $ smimixj - a)), j = 1, 2, • • • , M - 1, 



i=i 



n+l 



Si 



5 2 \r j *\ 2 + {5 1 -5 2 \r j *\ 2 )e- k ^ 



cxp 



-^ln 

2<5 2 



(2.44) 



Remark 2.4 // £/ie homogeneous peridic boundary conditions in (2.2) are replaced 
by the periodic boundary conditions 



i/>(a,t) =i/>(b,t), ip x (a,t) =4> x (b,t), t>0, 



(2.45) 



the TSSP scheme (2.13) still works provided that one replaces the sine-series in 
(2.13) by a Fourier-series [7, 6, 8]. 



2.3 Stability and decay rate 

Let U = {U^U^'-^UmY with U = U M 
/ 2 -norm on the interval (a, b), i.e., 



and 



be the usual discrete 



iifii. 



\ 



M \ 3\ 



(2.46) 



For the stability of the time-splitting sine-spectral approximations TSSP (2.13), 
we have the following lemma, which shows that the total normalization does not 
increase. 



Lemma 2.1 The time- splitting sine-spectral schemes (TSSP) (2.13) are uncondi- 
tionally stable if g(s) > for s > 0. In fact, for every mesh size h > and time 
step k > 0, 

U n+1 \\p < < = Uo\\p, n = 0, 1, 2, • • • (2.47) 

Furthermore, when a linear damping term is used in (1.5), i.e., we choose g{p) = 8 
with 5 > 0, the decay rate of the normalization satisfies 

= e- 2St "U% 2 = e - 2<5t "||^ || i2 , n = l,2,-.. (2.48) 

In fact, (2.48) is a discretized version of the decay rate of the normalization N(t) 
in (1.8). 
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Proof: We combine (2.13), (2.14), (2.46) and note that F(s,r) > for s > and 
r > 0, to obtain 



M-1 



b — a 



M 



M-1 



= m E^p[-2^(i^ri 2 ^/2)] 



1 M-1 

— y 

1 M-1 

- E 

2 ^ 

M-1 



M-1 



E 



i M-1 

< — E k 

1 M-1 



z=i 

> M-1 



2 1 ^ , 2 



/=i 



— E *HM X 3 - a)) 



3=1 



1 M-1 2 

i s k 



J'=l 



2 1 ._ 2 



- £ exp [-2F (|^| 2 , fe/2)] |^"| < - E \W 



M 



3=1 



1 



nil 2 



U 2 - 



b — a 

Here, we used the identity 

M-1 



(2.49) 



E sin 

3=1 



irr j 



sin 



M 



0, r - s ^ 2mM, 
= \M/2, r- S = 2mM,r^2nM, inte S er - 

(2.50) 

When a linear damping term is added to the NLS (1.5), the equality (2.48) 
follows from the above proof, Eq. (2.18), and 



M-1 

E ex p 

3=1 



M-1 



-2F(|^| 2 ,fc/2)] \^\ = E e- 5fc |^ n 



M-1 



3 Numerical examples 

In this section we present numerical tests of the TSSP (2.13) for solving a focusing 
cubic NLS appearing in nonlinear optics [18, 38] and for the Gross-Pitaeskii equation 
in BEC [8] in 2d with a linear, a cubic, or a quintic damping term. In our compu- 
tations, the initial condition (1.2) is always chosen such that |-^o( x )| decays to zero 
sufficiently fast as |x| — > oo. We choose an appropriately large rectangle [a, b] x [c, d] 
in 2d to avoid that the homogeneous periodic boundary condition (2.2) introduce 
a significant (aliasing) error relative to the whole space problem. To quantify the 
numerical results of the GPE for a BEC, we define the condensate widths along the 
x, y and z-axis by 



N(t) JR d 



J ^ a 2 \ip(-x., i)\ 2 <ix, with a = x, y, or z. 
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Example 1 Solution of the 2d damped focusing cubic nonlinear Schrodinger 
equation. We choose d — 2, a — 1 and V(x, y) = in (1.5) and present computations 
for three different damping terms (5 > 0): 

/. A linear damping term, i.e. we choose g(p) = 5. 

II. A cubic damping term, i.e. we choose g(p) = 5(3 p. 

III. A quintic damping term, i.e. we choose g(p) = 5f3 2 p 2 . 

The initial condition (1.6) is taken to be 

i,{x, y, 0) = Mx, y) = 4= e-^+^ 2 )/ 2£ , (*, y) e M 2 . (3.1) 

We assume 7y = 2, e = 0.2 and (3 = 8 in (1.5) such that £(0) = -0.751582 < in 
(1.4). We solve the NLS on the square [—16, 16] 2 , i.e., a = c = —16 and b = d = 16 
with mesh size h — ^, time step A; = 0.0002 and homogeneous periodic boundary- 
conditions along the boundary of the square. We compare the effect of changing the 
damping parameter S in the three different cases I, II and III. 

Figure 1 shows the surface plot of the density \ip(x, y, t)\ 2 at time t = 1.25 
with 5 = 0.5; plots of the normalization, energy and central density \ip(0, 0, t)\ 2 as 
functions of time with 5 = 0.5, 0.3 and 5 = (no damping) for case I. Figure 2 
shows similar results for case II and Figure 3 for case III. Furthermore Figure 4 
shows contour plots of the density |?/>| 2 at different times for case III with 5 = 0.01. 

In the numerically computations, a blowup is detected either from the plot of 
the central density 1^(0, 0, t)\ 2 which at the blowup shows a very sharp spike with 
a peak value that increases when the mesh size h decreases, or from the plot of the 
energy E(t) which has a very sharp spike with negative values at the blowup. In 
fact, the method TSSP (2.13) aims to capture the solution of damped NLS without 
blowup, i.e. physical revelant solution. If one wants to capture the blowup rate of 
NLS, we refer to [29, 33]. 

From the numerical results we find the following conditions for arresting a blowup 
of the wave function with initial energy E(0) < 0. (1) For linear damping the blowup 
is arrested if the damping parameter is bigger than a certain threshold value which 
we find to be <5 t h ~ 0.461 by numerical experiments. As shown in Fig. lb blowup is 
arrested for 5 = 0.5 > 5^ while the wave function blows up for 5 < 5 t h as can be seen 
from Fig. lc&d where we have chosen 5 = 0.3 < 5 th and 5 = < S th , respectively. 
The time at which the blowup of the wave function happens, however, increases with 
increasing 8 (cf. Fig. lc&d). (2) For a cubic damping term with 5 > the blowup 
of the wave function is always arrested (cf. Fig. 2). (3) The above observation (2) 
also holds for a quintic damping term (cf. Fig. 3). 

For linear damping, we also test the dependence of the threshold value of the 
damping parameter 8 t ^ on (3 and the initial data. First we take 7^ = 2 and e = 0.2 
in (3.1). Table 1 shows the threshold values <5 t h for different (3 in (1.5), and E(0) 
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(3 = 8 (3 = 16 /3 = 32 /3 = 64 /? = 128 



E(0) -0.7516 -5.253 -14.256 -32.263 -68.275 
5 th 0.461 3.655 10.35 22.15 40.05 



Table 1: Dependence of 5 t h on (5 for j y = 2 and £ = 0.2 in (3.1). 





e = 0.8 


e = 0.4 £ = 0.2 


e = 0.1 


£ = 0.05 


E(0) 


-1.3133 


-2.6266 -5.2532 


-10.506 


-21.013 




0.895 


1.845 3.655 


7.25 


14.55 



Table 2: Dependence of 5 t h on e in (3.1) for (5 — 16 in (1.5) and 7 y = 2 in (3.1). 

represents the initial energy. Then we choose /3 — 16 in (1.5) and 7 y = 2 in (3.1). 
Table 2 displays the threshold values 8 t h for different values of e in (3.1). 
From Table 1 we find by a least square fitting, 

5 th = -0.6930£(0) or S th = 0.3872/? - 2.4627. 

Similarly from Table 2 we obtain 

S th = -0.6922£(0). 

Based on this observation, we conclude that the threshold value of the linear damping 
parameter 5 th depends linearly on the initial energy E(0). 

Example 2 Solution of the 2d damped GPE with focusing nonlinearity. We 
choose d — 2, a — 1 and V(x,y) = \{^i1x 2 + 7^t/ 2 ) to be a harmonic oscillator 
potential with 7^,7^ > in (1.5). Again, we present computations for the same 
three different damping terms in (1.5) as those we studied in Example 1. 

We take 7^ = 1 and 7^ = 4. The initial condition (1.6) is assumed to be 
the ground-state solution of (1.5) with g{p) = (i.e. undamped case) and (5 = 
—40. The cubic nonlinearity is ramped linearly from (3 = —40 (defocusing) to 
(3 = 50 (focusing) during the time interval [0,0.1] and afterwards kept constant. 
The absorption parameter was set to 8 — during the time interval [0,0.1] and 
increased to a positive value 5 > afterwards. 

We solve the GPE on the rectangle [-24, 24] x [-6, 6], i.e., for a = -24, b = 24, 
c = — 6 and d — 6 with mesh size h x = ^, h y = yfg, time step k = 0.0005 and 
homogeneous periodic boundary conditions along the boundary of the rectangle. 
Again, we compare the effect of changing the damping parameter S in the three 
different cases I, II and III. 

Figure 5 shows a surface plot of the density \ip(x, y, t) | 2 at times t — (ground- 
state solution) and t = 2.8 with 5 = 1.25; normalization, energy and central density 
1^(0, 0, t)\ 2 as functions of time with 5 = 1.25, 1.1 and (no damping) for case I. 
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Figure 6 shows similar results for case II and Figure 7 for case III. Furthermore 
Figure 8 shows contour plots of the density \tp\ 2 at different times for case III with 
5 = 0.15. 

From our numerical results we find that the observations (l)-(3) made for exam- 
ple 1 are still valid with the additional trapping potential. However, the value of 5 t h, 
depends on (3 (or initial energy E(0)) and we find S th w 1.185 for linear damping 
(cf. Fig. 5). 

3.1 Discussion 

In this subsection we discuss our numerical results in terms of physical properties 
of a BEC described by the GPE. We concentrate on those cases where a collapse of 
the wave function is arrested since this collapse leads to unphysical processes like 
the negative peaks in the energy E(t) shown in Figs. lc&d,5e&f. 

The general form of the time evolution in example 1 is similar for all three cases. 
Initially the cloud of atoms contracts due to the attractive interaction between the 
particles. This contraction is accompanied by an increase in the energy due to 
particle loss which is most efficient in regions of high particle density. These regions 
are characterized by a negative local energy density leading to an increase in energy 
for each particle lost there. After the central particle density has reached a maximum 
the cloud starts to expand due to the kinetic energy gained by the particles during 
the contraction. Particles are emitted from the cloud in burst like pulses which can 
be seen in Figs. 4,8. Such bursts have also been seen in BEC experiments [13]. 
The main differences between the three cases are the behavior of the energy and 
the number of particles as a function of time. In case I where we assumed a linear 
damping term the loss rate of particles from the condensate is independent of the 
shape of the condensate wave function. The energy decrease during the condensate 
expansion is determined by the loss of particles (cf. Fig. lb). In the cases of cubic 
and quintic damping the loss term only has a significant effect on the time evolution 
of the condensate during the contraction. When the condensate expands the density 
of particles is so low that the loss terms have only a very small effect and the energy 
E(t) and the number of particles N(t) remain almost constant (see Figs. 2c,3c&d). 

In example 2 we add an additional trap potential which confines the BEC and 
assume a realistic scenario (described above) to prepare the condensate in the trap 
(cf. experiments by Donley et al. [13]). We find that the initial process of turning on 
the attractive interactions between the particles leads to oscillations in the widths of 
the condensate [8] as can be seen from Figs. 5,6,7. However, neither the additional 
trap potential nor these oscillations significantly alter the behavior of the system 
compared to example 1 when the condensate is strongly contracted. Before and 
after this contraction some differences can be seen. By looking at Figs. 5,6 we find 
that the first minimum in a y due to the oscillations of the condensate causes and 
increase in the central density and in the energy. For cubic and quintic damping this 
is accompanied by an increased particle loss. However, an arrested collapse of the 
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wave function only happens when both a x and a y attain a minimum value due to the 
attractive interactions (cf. Fig. 5d and Fig. 6b). We also note that the frequency of 
the oscillations after an arrested collapse has happened is not significantly influenced 
by the damping terms. The amplitude of these oscillations is, however, strongly 
dependent on 8 and decreases with increasing 5. Finally, we want to mention that a 
series of contractions and expansion of the condensate is possible. In Fig. 7b we find 
three contractions of the condensate where only the first one reaches a sufficiently 
high particle density to lead to an increase in energy while the next two contractions 
show a rather smooth decrease in energy and particle number. For a smaller quintic 
damping term we obtain two contractions of the condensate which increase the 
energy (see Fig. 7c). 

4 Conclusions 

We extended the explicit unconditionally stable second-order time-splitting sine- 
spectral (TSSP) method for solving damped focusing nonlinear Schrodinger equa- 
tions. We showed that this method is time transversal invariant and preserves the 
exact decay rate of the normalization for a linear damping of the NLS. Extensive 
numerical tests were presented for the cubic focusing nonlinear Schrodinger equation 
in 2d with linear, cubic and quintic damping terms. Our numerical results show that 
quintic damping always arrests blowup, whereas linear and cubic damping can arrest 
blowup only when the damping parameter 5 is bigger than a certain threshold value 
<5 t h- We will apply this novel method to solve the 3d Gross-Pitaevskii equation with 
a quintic damping term and compare the numerical results with the experimental 
dynamics [13] of collapsing and exploding BECs [9]. 
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Figure 1: Numerical results in Example 1 case I. a). Surface plot of the density 
|^| 2 at time t = 1.25 with 5 = 0.5. Normalization, energy and central density 
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Figure 2: Numerical results in Example 1 case II. Surface plot of the density \ip\ 2 
with 5 = 0.02: a). At time t = 0.4, b). t = 1.0. Normalization, energy and 
central density \ip(0,0,t)\ 2 as functions of time: c). with 5 = 0.02, d). 5 = 0.005 
(with h = 1/128, k = 0.00002). 
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c). t d). t 

Figure 3: Numerical results in Example 1 case III. Surface plot of the density \^\ 2 
with 5 = 0.01: a). At time t = 0.4, b). t = 1.0. Normalization, energy and central 
density |^(0,0,t)| 2 as functions of time: c). with 5 = 0.01, d). 5 = 0.001. 
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Figure 4: Contour plots of the density \ip\ 2 at different times in Example 1 case III 
with 5 = 0.01. a), t = 0, b). t = 0.2, c). * = 0.4, d). t = 0.6, e). t = 0.8, f). t = 1. 
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Figure 5: Numerical results in Example 2 case I. Surface plot of the density 
|-^| 2 with 5 = 1.25: a). At time t — (ground-state solution), b). t = 2.8. 
Normalization, energy and central density |?/>(0, 0, t)\ 2 as functions of time: c). with 
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Figure 6: Numerical results in Example 2 case II. a). Surface plot of the density 

\ip\ 2 with S = 0.15: At time t = 0.8 (left column) and t = 2.4 (right column). 

Normalization, energy and central density \ip(0, 0, t)| 2 (left column) and condensate 

widths (right column) as functions of time: b). With 5 = 0.15; c). 5 = 0.04 (with . 
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Figure 7: Numerical results in Example 2 case III. a). Surface plot of the density 
\ijj\ 2 with S = 0.15: At time t = 0.8 (left column) and t = 3.2 (right column). 
Normalization, energy and central density |^(0, 0, t)| 2 (left column) and condensate 
widths (right column) as functions of time: b). With 5 = 0.15; c). 5 = 0.005. 
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